#include <float.h>
#include <math.h>
#include <memory.h>
#include <stdio.h>
#include <stdlib.h>

#define KERNEL_LINEAR    //线性核函数
//#define KERNEL_RBF		//径向基函数核函数
//#define KERNEL_POLY	//多项式核函数

#define MAX_SIZE 8192
#define DIMISION	4
#define C			1

#define MAX_ITERATOR 0x1000


#define OFFSET(x, y) 	((x) > (y) ? (((x)+1)*(x) >> 1) + (y) : (((y)+1)*(y) >> 1) + (x))

typedef enum { POS = 1, NEG = -1 } label;

typedef struct SVM
{
	double weight[DIMISION];
	double bias;
}SVM;

typedef struct Dataset
{
	double data[MAX_SIZE][DIMISION];
	label label[MAX_SIZE];
	unsigned count;
}SVMDataset;

typedef struct
{
	SVM svm;
	SVMDataset dataset;
	double alpha[MAX_SIZE];
	double gradient[MAX_SIZE];
	double *cache;
}SVMBundle;

#ifdef KERNEL_LINEAR

double kernel(double src0[], double src1[], unsigned count)
{
	double result = 0;
	for (unsigned i = 0; i < count; ++i)
		result += src0[i] * src1[i];
	return result;
}

static SVMPrint(SVM *svm)
{
	printf("f(x)=%lf", svm->bias);
	for (unsigned i = 0; i < DIMISION; ++i)
		printf("+%lf*x%d", svm->weight[i], i);
	printf("\n");
}

#else
#ifdef KERNEL_RBF

#define GAMMA 0.1
double kernel(double src0[], double src1[], unsigned count)
{
	double result = 0;
	for (unsigned i = 0; i < count; ++i)
		result += (src0[i] - src1[i]) * (src0[i] - src1[i]);
	return exp(-result * GAMMA);
}

static SVMPrint(SVM *svm)
{
	printf("f(x)=%lf+exp(%lf*(", svm->bias, GAMMA);
	for (unsigned i = 0; i < DIMISION; ++i)
	{
		printf("(%lf-x%d)^2", svm->weight[i], i);
		if (i < DIMISION - 1)
			printf(" + ");
		else
			printf("))\n");
	}
}
#undef GAMMA

#else
#ifdef KERNEL_POLY

#define K 2
double kernel(double src0[], double src1[], unsigned count)
{
	double result = 0;
	for (unsigned i = 0; i < count; ++i)
	{
		double x0 = 1, x1 = 1;
		for (unsigned j = 0; j < K; ++j)
		{
			x0 *= src0[i];
			x1 *= src1[i];
			result += x0 * x1;
		}
	}
	return result;
}

static SVMPrint(SVM *svm)
{
	printf("f(x) = %lf", svm->bias);
	for (unsigned i = 0; i < DIMISION; ++i)
	{
		double weight = svm->weight[i];
		printf("+%lf*x%d", weight, i);
		for (unsigned j = 1; j < K; ++j)
		{
			weight *= svm->weight[i];
			printf("+%le*x%d^%d", weight, i, j + 1);
		}
	}
	printf("\n");
}
#undef K

#else
#error "Kernel Not Define!"
#endif // KERNEL_POLY
#endif // KERNEL_RBF
#endif // KERNEL_LINEAR


void SVMInitial(SVMBundle *bundle)
{
	bundle->cache = (double *)calloc((MAX_SIZE + 1)*MAX_SIZE >> 1, sizeof(double));
	double *cache = bundle->cache;
	double *alpha = bundle->alpha;
	double *gradient = bundle->gradient;
	double (*x)[DIMISION] = bundle->dataset.data;
	label *y = bundle->dataset.label;
	unsigned count = bundle->dataset.count;
	for (unsigned i = 0; i < count; ++i)
		for (unsigned j = 0; j <= i; ++j)
			cache[OFFSET(i, j)] = y[i] * y[j] * kernel(x[i], x[j], DIMISION);
	for (unsigned i = 0; i < count; ++i)
	{
		gradient[i] = -1;
		for (unsigned j = 0; j < count; ++j)
			gradient[i] += cache[OFFSET(i, j)] * alpha[j];
	}
}

void SVMDispose(SVMBundle *bundle)
{
	free(bundle->cache);
	bundle->cache = 0;
}

//static int loadTrainset(SVMBundle *bundle, double dataset[DIMISION], label result)
//{
//	unsigned count = bundle->dataset.count, index = count;
//	double(*data)[DIMISION] = bundle->dataset.data;
//	label *label = bundle->dataset.label;
//	if (count < MAX_SIZE)
//	{
//		count = ++bundle->dataset.count;
//	}
//	else
//	{
//		double *weight = bundle->svm.weight;
//		double bias = bundle->svm.bias;
//		double maxmargin = kernel(weight, dataset, DIMISION) + bias;
//		if (maxmargin < 0) maxmargin = -maxmargin;
//		for (unsigned i = 0; i < MAX_SIZE; ++i)
//		{
//			double margin = kernel(weight, data[i], DIMISION) + bias;
//			if (margin < 0) margin = -margin;
//			if (margin > maxmargin)
//			{
//				maxmargin = margin;
//				index = i;
//			}
//		}
//		if (index == count) return 0;
//	}
//	memcpy(data[index], dataset, sizeof(data[index]));
//	label[index] = result;
//	for (unsigned i = 0; i < count; ++i)
//		bundle->cache[OFFSET(index, i)] = label[i] * result * kernel(dataset, data[i], DIMISION);
//	return 1;
//}


static int SVMSolve(SVMBundle *bundle)
{
	unsigned count = bundle->dataset.count;
	double *gradient = bundle->gradient;
	double *cache = bundle->cache;
	double *alpha = bundle->alpha;
	double *weight = bundle->svm.weight;
	double(*data)[DIMISION] = bundle->dataset.data;
	label *y = bundle->dataset.label;

	{
		unsigned posCount = 0;
		for (unsigned i = 0; i < count; ++i)
			posCount += y[i] == POS;
		if (posCount == 0 || posCount == count)
			return 1;
	}

	for (unsigned i = 0; i < MAX_ITERATOR; ++i)
	{
		unsigned x0 = 0, x1 = 1;
		//根据梯度选取进行优化的alpha值
		{
			double gmax = -DBL_MAX, gmin = DBL_MAX;
			for (unsigned i = 0; i < count; ++i)
			{
				if ((alpha[i] < C && y[i] == POS || alpha[i] > 0 && y[i] == NEG) && -y[i] * gradient[i] > gmax)
				{
					gmax = -y[i] * gradient[i];
					x0 = i;
				}
				else if ((alpha[i] < C && y[i] == NEG || alpha[i] > 0 && y[i] == POS) && -y[i] * gradient[i] < gmin)
				{
					gmin = -y[i] * gradient[i];
					x1 = i;
				}
			}
		}

		double oldAlpha0 = alpha[x0];
		double oldAlpha1 = alpha[x1];
		//调整选取的alpha值
		if (y[x0] != y[x1])
		{
			double coef = cache[OFFSET(x0, x0)] + cache[OFFSET(x1, x1)] + 2 * cache[OFFSET(x0, x1)];
			if (coef <= 0) coef = DBL_MIN;
			double delta = (- gradient[x0] - gradient[x1]) / coef;
			double diff = alpha[x0] - alpha[x1];
			alpha[x0] += delta;
			alpha[x1] += delta;
			unsigned max = x0, min = x1;
			if (diff < 0)
			{
				max = x1;
				min = x0;
				diff = -diff;
			}
			if (alpha[max] > C)
			{
				alpha[max] = C;
				alpha[min] = C - diff;
			}
			if (alpha[min] < 0)
			{
				alpha[min] = 0;
				alpha[max] = diff;
			}
		}
		else
		{
			double coef = cache[OFFSET(x0, x0)] + cache[OFFSET(x1, x1)] - 2 * cache[OFFSET(x0, x1)];
			if (coef <= 0) coef = DBL_MIN;
			double delta = (-gradient[x0] + gradient[x1]) / coef;
			double sum = alpha[x0] + alpha[x1];
			alpha[x0] += delta;
			alpha[x1] -= delta;
			unsigned max = x0, min = x1;
			if (alpha[x0] < alpha[x1])
			{
				max = x1;
				min = x0;
			}
			if (alpha[max] > C)
			{
				alpha[max] = C;
				alpha[min] = sum - C;
			}
			if (alpha[min] < 0)
			{
				alpha[min] = 0;
				alpha[max] = sum;
			}
		}
		double delta0 = alpha[x0] - oldAlpha0;
		double delta1 = alpha[x1] - oldAlpha1;
		if (DBL_MIN > (delta0 < 0 ? -delta0 : delta0) && DBL_MIN >(delta1 < 0 ? -delta1 : delta1))
			break;
		for (unsigned i = 0; i < count; ++i)
			gradient[i] += cache[OFFSET(i, x0)] * delta0 + cache[OFFSET(i, x1)] * delta1;
	}

	//计算权重
	{
		memset(weight, 0, sizeof(double)*DIMISION);
		for (unsigned j = 0; j < DIMISION; ++j)
			for (unsigned i = 0; i < count; ++i)
				if (alpha[i] != 0)
					weight[j] += alpha[i] * y[i] * data[i][j];
	}
	//计算偏移常量
	{
		double maxneg = -DBL_MAX, minpos = DBL_MAX;
		SVM *svm = &bundle->svm;
		for (unsigned i = 0; i < count; ++i)
		{
			double wx = kernel(svm->weight, data[i], DIMISION);
			if (y[i] == POS && minpos > wx)
				minpos = wx;
			else if (y[i] == NEG && maxneg < wx)
				maxneg = wx;
		}
		svm->bias = -(minpos + maxneg) / 2;
	}
	return 0;
}



int LoadIrisDataset(SVMDataset *dataset)
{
	FILE *fp = fopen("iris.csv", "r");
	if (!fp) return 1;
	unsigned *count = &dataset->count;
	char temp[32];
	fscanf(fp, "%u,%s\n", count, temp);
	for (unsigned i = 0; i < *count; ++i)
	{
		double *data = dataset->data[i];
		int type = 0;
		fscanf(fp, "%lf,%lf,%lf,%lf,%d\n", data, data + 1, data + 2, data + 3, &type);
		dataset->label[i] = type > 0 ? POS : NEG;
	}
	fclose(fp);
	return 0;
}




int main()
{
	SVMBundle bundle = { 0 };
	SVMDataset *dataset = &bundle.dataset;
	if (LoadIrisDataset(&bundle.dataset))
	{
		printf("Dataset iris.csv Not Found!\n");
		return 1;
	}
	SVMInitial(&bundle);
	SVMSolve(&bundle);
	SVMDispose(&bundle);

	unsigned right = 0;
	SVM *svm = &bundle.svm;
	for (unsigned i = 0; i < dataset->count; ++i)
	{
		double margin = kernel(svm->weight, dataset->data[i], DIMISION) + svm->bias;
		if (dataset->label[i] == POS)
		{
			right += margin > 0;
			//printf("POS\tMargin:%le\n", margin);
		}
		else
		{
			right += margin < 0;
			//printf("NEG\tMargin:%le\n", margin);
		}
	}
	printf("Right:%u/%u\n", right, dataset->count);
	SVMPrint(&bundle.svm);
	system("pause");
}